//////////////////////////////////////////////////////////////////
//                        COPYRIGHT NOTICE
// (c) Copyright 2003- Jeongnim Kim
//   All rights reserved.
//////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////
//                    CONTACT INFORMATION
//   Jeongnim Kim
//   National Center for Supercomputing Applications &
//   Materials Computation Center
//   University of Illinois, Urbana-Champaign
//   Urbana, IL 61801
//   e-mail: jnkim@ncsa.uiuc.edu
//   Tel:    217-244-6319 (NCSA) 217-333-3324 (MCC)
//
// Supported by 
//   National Center for Supercomputing Applications, UIUC
//   Materials Computation Center, UIUC
//   Department of Physics, Ohio State University
//   Ohio Supercomputer Center
//////////////////////////////////////////////////////////////////
// -*- C++ -*-

#ifndef OHMMS_GENERIC_PARPACKSYMENGINE_H
#define OHMMS_GENERIC_PARPACKSYMENGINE_H

#include <complex>
#include <iostream>
#include <map>
#include <algorithm>
using namespace std;
#include "Message/CommOperators.h"
#include "EigSolver/EigenSpaceDimension.h"
#include "OhmmsNumeric/OhmmsBlasf.h"
#include "OhmmsPETE/OhmmsVector.h"
#ifdef ADD_
#define pssaupd pssaupd_
#define psseupd psseupd_
#define pdsaupd pdsaupd_
#define pdseupd pdseupd_
#define pznaupd pznaupd_
#define pzneupd pzneupd_
#endif
extern "C" {

  void pssaupd(int *comm,
	       int *ido, const char *bmat, int *n, const char *which,
	       int *nev, float *tol, float *resid,
	       int *ncv, float *V, int *ldv,
	       int *iparam, int *ipntr, float *workd,
	       float *workl, int *lworkl, int *info);
  
  void pdsaupd(int *comm,
	       int *ido, const char *bmat, int *n, const char *which,
	       int *nev, double *tol, double *resid,
	       int *ncv, double *V, int *ldv,
	       int *iparam, int *ipntr, double *workd,
	       double *workl, int *lworkl, int *info);

  void pznaupd(int *comm,
	       int *ido, const char* bmat, int *n, const char* which,
	       int *nev, double *tol, complex<double> *resid,
	       int *ncv, complex<double> *V, int *ldv,
	       int *iparam, int *ipntr, complex<double> *workd,
	       complex<double> *workl, int *lworkl, double *rwork,
	       int *info); 

  void psseupd(int *comm,
	       int *rvec, char *HowMny, int *select,
	       float *d, float *Z, int *ldz,
	       float *sigma, const char *bmat, int *n,
	       const char *which, int *nev, float *tol,
	       float *resid, int *ncv, float *V,
	       int *ldv, int *iparam, int *ipntr,
	       float *workd, float *workl,
	       int *lworkl, int *info);

  void pdseupd(int *comm,
	       int *rvec, char *HowMny, int *select,
	       double *d, double *Z, int *ldz,
	       double *sigma,  const char *bmat, int *n,
	       const char *which, int *nev, double *tol,
	       double *resid, int *ncv, double *V,
	       int *ldv, int *iparam, int *ipntr,
	       double *workd, double *workl,
	       int *lworkl, int *info);
  
  void pzneupd(int *comm,
	       int *rvec, const char* HowMny, int *select,
	       complex<double> *d, complex<double> *Z, int *ldz,
	       complex<double> *sigma, complex<double>* workev, 
	       const char* bmat, int *n,
	       const char* which, int *nev, double *tol, 
	       complex<double>* resid, int* ncv, complex<double> *V,
	       int *ldv, int *iparam, int *ipntr,
	       complex<double> *workd, complex<double> * workl,
	       int *lworkl, double *rwork, int *info); 
  
}

template<class T, class EigVec_t>
struct ParpackSymSolver {

  typedef typename OHMMS::type_traits<T>::value_type value_type;
  typedef typename OHMMS::type_traits<T>::real_type  real_type;
  typedef EigVec_t                 EigenStates_t;
  typedef Vector<value_type> WorkArray_t;
  typedef Vector<real_type>  RealArray_t;
  typedef Vector<int>        IndexArray_t;

  typedef typename WorkArray_t::size_type   size_type;

  int N, Nloc;
  int M;
  int ncv; // number of Lanczos vectors
  int info;
  int lworkl;
  int ido;
  int ierr;
  int nconv;
  int maxiter;
  int ishfts;
  int mode;
  int rvec;
  int comm;
  string which;
  int iparam[11], ipntr[14];
  char bmat;

  IndexArray_t select;
  IndexArray_t rcounts, displs;
  WorkArray_t gV; // global array to perform apply
  WorkArray_t workd; // array of 3N+1
  WorkArray_t workl; // array of lworkl
  WorkArray_t resid; // residue
  WorkArray_t eigVal;// eigen values
  WorkArray_t eigVec;// eigen vectors 
  RealArray_t rwork; // only for a complex problem

  real_type tol;  // tolerance
  value_type sigma; // reference value for shift-invert


  ParpackSymSolver(EigenSpaceDimension& dim, bool nonortho=false){   
    N = 0;
    M = 0;
    tol = 1.0e-13;
    info = 0;
    ido = 0;
    maxiter = 1000;
    ncv = -1;
    iparam[0] = 1; iparam[6] = 1;
    sigma = 0.0;
    which = "SM";
    ///get the communicator
    comm = Controller->getID();
  }

  ~ParpackSymSolver() { }

  void gatherv(double* pl, double *pg) {
    MPI_Allgatherv(pl, Nloc,MPI_DOUBLE, pg, rcounts.data(), displs.data(), 
		   MPI_DOUBLE,comm);
  }

  template<class HamiltonianMatrix_t>
  bool solve(HamiltonianMatrix_t& H) {

    bmat = 'I';
    Nloc =H.numcol(0);
    initializeWorkSpace();
    iparam[0] = 1; iparam[6] = 1;
    ido = 0;

    do{
      solve(eigVec.data());
      gV = 0.0;
      gatherv(workd.data()+ipntr[0]-1,gV.data());
      H.apply(gV.data(),workd.data()+ipntr[1]-1,true);
      //H.apply(workd.data()+ipntr[0]-1,workd.data()+ipntr[1]-1,false);
    } while(ido == -1 || ido == 1);

    if(info < 0) {
      cerr << endl << "Error with __aupd, info = " <<  info << endl;
      cerr << " Check documentation in __aupd " << endl << endl;
      return false;
    } else {
      rvec = 1;
      finalize(eigVec.data());
      return true;
    }
  }

  template<class HamiltonianMatrix_t>
  bool solve(HamiltonianMatrix_t& H, EigenStates_t& Psi) {
    bool done = solve(H);
    // copy eigen states according to the ascending eigenvalues
    multimap<double,int> ordered_e;
    for(int i=0; i<M; i++) ordered_e.insert(pair<double,int>(eigVal[i],i));
    multimap<double,int>::iterator it = ordered_e.begin();
    int ieig = 0;
    for(; it!=ordered_e.end(); it++) {
      eigVal[ieig] = (*it).first;      
      gatherv(eigVec.data()+(*it).second*Nloc,Psi.data()+ieig*N);
      //int is = (*it).second*Nloc;
      //std::copy(eigVec.data()+is,eigVec.data()+is+N, Psi.data()+ieig*N);
      ieig++; 
    }
    return done;
  }

  template<class HamiltonianMatrix_t>
  bool solve(HamiltonianMatrix_t& H, const HamiltonianMatrix_t& S) {
    return false;
  }

  template<class HamiltonianMatrix_t>
  bool solve(HamiltonianMatrix_t& H, const HamiltonianMatrix_t& S, 
	     EigenStates_t& Psi) {
    return false;
  }

  template<class HamiltonianMatrix_t>
  bool checkEigvectors(HamiltonianMatrix_t& H, ostream& os) {

    if ( ierr != 0) {
      os << endl << "Error with __eupd, info = " << ierr << endl;
      os << "Check the documentation of __eupd. " << endl;
      return false;
    } else {
      
      nconv =  iparam[4];
      if(rvec == 1) {
	WorkArray_t ax(Nloc);
	WorkArray_t res(M);
	
	for(int j=0; j<nconv; j++) { // residual calculations
	  ax = 0.0;
          gV = 0.0;
          gatherv(eigVec.data()+j*Nloc,gV.data());
	  H.apply(gV.data(), ax.data(),true);
	  value_type da  = -eigVal[j];
	  int incx = 1;
	  int incy = 1;
	  OhmmsBlas::axpy(N, da, eigVec.data()+j*Nloc, incx, ax.data(), incy);  
	  res[j] = OhmmsBlas::dot(Nloc, ax.data(), ax.data());
          gsum(res[j],comm);
	  //res[j] = OhmmsBlas::norm2(N, ax.data(), incx);
	  res[j] = sqrt(res[j])/fabs(eigVal[j]);
	}
	os << "Ritz values and relative residuals " << endl;
	for(int j=0; j<nconv; j++)
	  os << j << "  "  << eigVal[j] << "  " << res[j] << endl;
      } else {
	os << "Ritz values only " << endl;
	for(int j=0; j<nconv; j++)
	  os << j << "  "  << eigVal[j] <<  endl;
      }

      os << endl;
      if(info == 1) {
	os << endl << " Maximum number of iterations reached." << endl;
      } else if ( info == 3) {
	os << endl << "No shifts could be applied during implicit "
	   <<   "  Arnoldi update, try increasing NCV." << endl;
      }
      os << endl;
      os << " _SDRV1 " << endl;
      os << " ====== " << endl << endl;
      os << " Size of the matrix is " << N << endl;
      os << " The number of Ritz values requested is " << M << endl;
      os << " The number of Arnoldi vectors generated  (NCV) is " << ncv << endl;
      os << " What portion of the spectrum: " << which << endl;
      os << " The number of converged Ritz values is " << nconv << endl;
      
      os << " The number of Implicit Arnoldi update iterations taken is ";
      os <<  iparam[2] << endl;
      os << " The number of OP*x is " << iparam[8] << endl;
      os << " The convergence criterion is " << tol << endl;
      return true;
    }
  }

  template<class HamiltonianMatrix_t>
  bool checkEigvectors(const HamiltonianMatrix_t& H, 
		       const HamiltonianMatrix_t& S,
		       ostream& os) {
    return false;
  }

  inline void setMaxIter(int n) { maxiter = n;}
  inline void setMaxTrialStates(int n) { ncv = n;}
  inline void setTolerance(real_type x) { tol = x;}
  inline real_type eigval(size_type i) const { return eigVal[i];}

  //bool put(xmlDocPtr doc, xmlNsPtr ns, xmlNodePtr cur);
  bool put(xmlDocPtr doc,xmlNodeSetPtr nodeset);

  void initializeWorkSpace() {

    int node = Controller->ncontexts();
    int myid = Controller->mycontext();
    rcounts.resize(node);
    displs.resize(node+1);
    MPI_Allgather(&Nloc, 1, MPI_INT, rcounts.data(), 1, MPI_INT, comm);
    displs[0] = 0;
    for(int i=0; i<node; i++) displs[i+1] = displs[i]+rcounts[i];


    if(ncv < 0) ncv = 4*M;
    if(ncv > N) ncv = N;

    //TODO: COMPLEX and DOUBLE  (maybe not)
    lworkl = ncv*(ncv+8);
    iparam[2] = maxiter;
    if(workd.size() < 3*Nloc) {
      workd.resize(3*Nloc);
    }                                    
    if(resid.size() < Nloc) {
      resid.resize(Nloc);
    }
    if(select.size() < ncv) {
      select.resize(ncv);
    }
    if(eigVal.size() < ncv*2) {
      eigVal.resize(ncv*2);
    }
    if(eigVec.size() < ncv*Nloc) {
      eigVec.resize(ncv*Nloc);
    }
    if(workl.size() < lworkl) {
      workl.resize(lworkl);
    }
    if(gV.size()<N) {
      gV.resize(N);
    }    
  }

  void solve(double* eigv){
    int ldv = Nloc; 
    int nev = M;
    pdsaupd(&comm,
	    &ido, &bmat, &Nloc, which.c_str(), &nev, &tol,
	    resid.data(), &ncv, eigv, &ldv, iparam, ipntr,
	    workd.data(), workl.data(), &lworkl, &info);
  }

  void solve(complex<double>* eigv){

    if(rwork.empty()) {
      rwork = RArray_t(ncv); // only for a complex problem
    } 
    int ldv = Nloc; 
    int nev = M;
    pznaupd(&comm,
	    &ido, &bmat, &Nloc, which.c_str(), &nev, &tol,
	    resid.data(), &ncv, eigv, &ldv, iparam, ipntr,
	    workd.data(), workl.data(), &lworkl, rwork.data(), 
	    &info);
  }

  void solve(float* eigv){
    int ldv = Nloc; 
    int nev = M;
    pssaupd(&comm,
	    &ido, &bmat, &Nloc, which.c_str(), &nev, &tol,
	    resid.data(), &ncv, eigv, &ldv, iparam, ipntr,
	    workd.data(), workl.data(), &lworkl, &info);
  }

  void finalize(double *eigv) {
    char howmany ='A';
    int ldv = Nloc; 
    int nev = M;
    pdseupd( &comm,
	     &rvec, &howmany, select.data(),
	     eigVal.data(), eigv, &ldv, &sigma,
	     &bmat, &Nloc, which.c_str(), &nev, &tol, resid.data(), &ncv,
	     eigv, &ldv, iparam, ipntr,
	     workd.data(), workl.data(), &lworkl, &ierr );
    
    //for(int i=0; i<States.M; i++) ene[i] = eigVal[i];
  }

  void finalize(complex<double> *eigv) {
    char howmany ='A';
    int ldv = Nloc; 
    int nev = M;
    WorkArray_t workev(2*ncv);    
    pzneupd(&comm,
	    &rvec, &howmany, select.data(),
            eigVal.data(), eigv, &ldv, &sigma, 
            workev.data(),  &bmat, &Nloc, which.c_str(), &nev, &tol, 
            resid.data(), &ncv, eigv, &ldv, iparam, ipntr, workd.data(), 
            workl.data(), &lworkl, rwork.data(), &ierr );
  }

  void finalize(float *eigv) {
    char howmany ='A';
    int ldv = Nloc; 
    int nev = M;
    psseupd(&comm,
	    &rvec, &howmany, select.data(),
            eigVal.data(), eigv, &ldv, &sigma,
            &bmat, &Nloc, which.c_str(), &nev, &tol, resid.data(), &ncv,
            eigv, &ldv, iparam, ipntr,
            workd.data(), workl.data(), &lworkl, &ierr );
  }
};
#endif

/***************************************************************************
 * $RCSfile: ParpackSymSolver.h,v $   $Author: jnkim $
 * $Revision: 1.1.1.1 $   $Date: 2004/07/23 21:04:02 $
 * $Id: ParpackSymSolver.h,v 1.1.1.1 2004/07/23 21:04:02 jnkim Exp $ 
 ***************************************************************************/

